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Abstract. In horizontally shaken granular material different types of pattern formation have been reported. 
We want to deal with the convection instability which has been observed in experiments and which recently 
has been investigated numerically. Using two dimensional molecular dynamics we show that the convection 
pattern depends crucially on the inelastic properties of the material. The concept of restitution coefficient 
provides arguments for the change of the behaviour with varying inelasticity. 
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1 Introduction 



When granular material in a cubic container is shaken 
horizontally one observes experimentally different types of 
instabilities, i.e. spontaneous formation of ripples in shal- 
low beds jjj, liquefaction 0,0], convective motion [|[|^] and 
recurrent swelling of shaken material where the period of 
swelling decouples from the forcing period |^|. Other in- 
teresting experimental results concerning simultaneously 
vertically and horizontally vibrated granular systems []7J 
and enhanced packing of spheres due to horizontal vibra- 
tions II have been reported recently. Horizontally shaken 
granular systems have been simulated numerically using 
cellular automata M as well as molecular dynamics tech- 
niques |,§,|M|l!Q. Theoretical work on horizontal 
shaking can be found in JTT| and the dynamics of a single 
particle in a horizontally shaken box has been discussed 
in H. 




Recently the effect of convection in a horizontally shaken 
box filled with granular material attracted much atten- 
tion and presently the effect is studied experimentally by 
different groups 0,0,^). Unlike the effect of convective 
motion in vertically shaken granular material which has 
been studied intensively experimentally, analytically and 
by means of computer simulations (s. e.g. ]l4] , |l^ ,|l^,|l7|), 
there exist only a few references on horizontal shaking. 
Different from the vertical case, where the "architecture" 
of the convection pattern is very simple |lS| , in horizon- 
tally shaken containers one observes a variety of different 
patterns, convecting in different directions, in parallel as 
well as perpendicular to the direction of forcing Q . Under 
certain conditions one observes several convection rolls on 
top of each other ||. An impression of the complicated 



convection can be found in the internet |19 



Fig. 1. Sketch of the simulated system. 



Whereas the properties of convection in vertically sha- 
ken systems can be reproduced by two dimensional molec- 
ular dynamics simulations with good reliability, for the 
case of horizontal motion the results of simulations are in- 
consistent with the experimental results: in all experimen- 
tal investigations it was reported that the material flows 
downwards close to the vertical walls Q,|^,^,^9|, but re- 
ported numerical simulations systematically show surface 
rolls in opposite direction accompanying the more realistic 
deeper rolls, or even replacing them completely |jic|| . 

Our investigation is thus concerned with the convec- 
tion pattern, i.e. the number and direction of the convec- 
tion rolls in a two dimensional molecular dynamics sim- 
ulation. We will show that the choice of the dissipative 
material parameters has crucial influence on the convec- 
tion pattern and, in particular, that the type of convection 
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rolls observed experimentally can be reproduced by using 
sufficiently high dissipation constants. 



2 Numerical Model 

The system under consideration is sketched in Fig. |l]: 
we simulate a two-dimensional vertical cross section of 
a three-dimensional container. This rectangular section 
of width L = 100 (all units in cgs system), and infinite 
height, contains N = 1000 spherical particles. The sys- 
tem is periodically driven by an external oscillator x(t) = 
A sin(27r/t) along a horizontal plane. For the effect we 
want to show, a working frequency / = 10 and amplitude 
A = 4 is selected. These values give an acceleration am- 
plitude of approximately 1617. Lower accelerations affect 
the intensity of the convection but do not change the ba- 
sic features of the convection pattern which we want to 
discuss. As has been shown in past the fluidization 
point, a much better indicator of the convective state is 
the dimensionless velocity A2nf / \fLg~. This means that in 
small containers motion saturates earlier, hence, results 
for different container lengths at the same values of the 
acceleration amplitude cannot be compared directly. Our 
acceleration amplitude ~ Vog corresponds to ~ 3g in a 
10 cm container (provided that the frequency is the same 
and particle sizes have been scaled by the same amount). 

The radii of the particles of density 2 are homoge- 
neously distributed in the interval [0.6,1.4]. The rough 
inner walls of the container are simulated by attaching ad- 
ditional particles of the same radii and material properties 
(this simulation technique is similar to "real" experiments, 
e.g. H). 

For the molecular dynamics simulations, we apply a 
modified soft-particle model by Cundall and Strack fl20f : 
Two particles i and j, with radii Ri and Rj and at posi- 
tions Yi and Tj, interact if their compression £y = Rj, + 
Rj — |r s ; — Vj \ is positive. In this case the colliding spheres 
feel the force n N + FfjYi s , with n" and n s being the 
unit vectors in normal and shear direction. The normal 
force acting between colliding spheres reads 



Hj = {miirij) j (rrii 
the relative velocity at the point of contact 




where Y is the Young modulus, v is the Poisson ratio and 
B is a material constant which characterizes the dissipa- 
tive character of the material [^l| . 

Rf=(R i R j )/(R i + R j ) (2) 

is the effective radius. For a strict derivation of ([j]) see ^l[ 

For the shear force we apply the model by Haff and 
Werner |23| 
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rrij ) and 



(4) 



F§ = sign (v[f^j min { 
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, rel 



(3) 



fli and flj are the angular velocities of the particles. 

The resulting momenta Mi and Mj acting upon the 
particles are Mi = F?Ri and Mj = -F$Rj. Eq. (|) takes 
into account that the particles slide upon each other for 
the case that the Coulomb condition fj, \ F^ \ < \F$\ 
holds, otherwise they feel some viscous friction. By means 
of 7„ = BY/(\ — v 2 ) and 7 S , normal and shear damping 
coefficients, energy loss during particle contact is taken 
into account p4[ . 

The equations of motion for translation and rotation 
have been solved using a Gear predictor-corrector scheme 
of sixth order (e.g. p5[). 

The values of the coefficients used in simulations are 
y/(l - v 2 ) = 1 X 10 8 , 7 S = 1 x 10 3 , ju = 0.5. For the effect 
we want to show, the coefficient 7„ takes values within the 
range [l0 2 ,10 4 ]. 



3 Results 

The mechanisms for convection under horizontal shak- 
ing have been discussed in pQ] . Now we can show that 
these mechanisms can be better understood by taking 
into account the particular role of dissipation in this prob- 
lem. The most striking consequence of varying the normal 
damping coefficient is the change in organization of the 
convective pattern, i.e. the direction and number of rolls 
in the stationary regime. This is shown in Fig. |^, which 
has been obtained after averaging particle displacements 
over 200 cycles (2 snapshots per cycle) . The asymmetry of 
compression and expansion of particles close to the walls 
(where the material results highly compressible) explains 
the large transverse velocities shown in the figure. Note, 
however, that the upward and downward motion at the 
walls cannot be altered by this particular averaging pro- 
cedure. 

The first frame shows a convection pattern with only 
two rolls, where the arrows indicate that the grains slide 
down the walls, with at most a slight expansion of the 
material at the surface. There are no surface rolls. This is 
very similar to what has been observed in experiments^,^, 
^| . In this case, dissipation is high enough to damp most of 
the sloshing induced by the vertical walls, and not even the 
grains just below the surface can overcome the pressure 
gradient directed downwards. 

For lower damping, we see the developing of surface 
rolls, which coexist with the inner rolls circulating in the 
opposite way. Some energy is now available for upward mo- 
tion when the walls compress the material fluidized during 
the opening of the wall "gap" (empty space which is cre- 
ated alternatively during the shaking motion) . This is the 
case reported in jnj . The last frames demonstrate how the 
original rolls vanish at the same time that the surface rolls 
grow occupying a significant part of the system. Another 
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feature shown in the figure is the thin layer of material in- 
volving 3 particle rows close to the bottom, which perform 
a different kind of motion. This effect, which can be seen 
in all frames, is due to the presence of the constraining 
boundaries but has not been analyzed separately. 
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Fig. 2. Velocity field obtained after cycle averaging of particle displacements, for different values of the normal damping 
coefficient, ■y n . The first one is 1 x 10 4 , and for obtaining each subsequent frame the coefficient has been divided by two. The 
frames are ordered from left to right and from top to bottom. The cell size for averaging is approximately one particle diameter. 
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Fig. 3. Horizontal (left) and vertical (right) displacements at 
selected positions of the frames in Fig. y (see the text for de- 
tails), for decreasing normal damping and as a function of 
depth. White indicates strongest flow along positive axis di- 
rections (up,right), and black the corresponding negative ones. 
The black region at the bottom of the left picture corresponds 
to the complex boundary effect observed in Fig. ^, involving 
only two particle layers. The figure below shows a typical con- 
vection pattern together with the sections at x = 10 and x = 30 
at which the displacements were recorded. 



With decreasing normal damping 7„ there are two 
transitions observable in Fig. ^, meaning that the convec- 
tion pattern changes qualitatively at these two particular 
values of j n : The first transition leads to the appearance 
of two surface rolls laying on top of the bulk cells and 
circulating in opposite direction. The second transition 
eliminates the bulk rolls. A more detailed analysis of the 
displacement fields (Fig. ||) allows us to locate the transi- 
tions much more precisely. In Fig. |§| we have represented in 
grey-scale the horizontal and vertical components of the 
displacement vectors pictured in Fig. ^| but in a denser 
sampling, analyzing data from 30 simulations correspond- 
ing to values of the normal damping coefficient within the 
interval [50,10000]. For horizontal displacements, we have 
chosen vertical sections at some representative position in 
horizontal direction (x — 30). For the vertical displace- 
ments, vertical sections of the leftmost part of the con- 
tainer were selected (x = 10), s. Fig. |j| lower part. 

The horizontal axis shows the values of the normal 
damping coefficient scaled logarithmically in decreasing 
sequence. The vertical axis represents the position in verti- 
cal direction, with the free surface of the system located at 



y w 60. One observes first that white surface shades, com- 
plemented by subsurface black ones, appear quite clearly 
at about 7 =2000 in Fig. || (left), indicating the appear- 
ance of surface rolls. On the other hand, Fig. || (right) 
shows a black area (indicative of downward flow along 
the vertical wall) that vanishes at 7„ ~ 200 (at this point 
the grey shade represents vanishing vertical velocity) . The 
dashed lines in Fig. || lead the eye to identify the transi- 
tion values. In the interval 200 < j n < 2000 surface and 
inner rolls coexist, rotating in opposite directions. 

One can analyze the situation in terms of the resti- 
tution coefficient. From Eq. (Q), the equation of motion 
for the displacement £y can be integrated and the relative 
energy loss in a collision rj = (Eq —E)/Eq (with E and Eq 
being the energy of the relative motion of the particles) 
can be evaluated approximately. Up to the lowest order in 
the expansion parameter, one finds [Era 



rj = 1.78 



1/5 



(5) 



where vq is the relative initial velocity in normal direc- 
tion, and r, I, time and length scales associated with the 
problem (see pqj for details), 



-B 




(G) 



For 7„ = 10 4 (the highest value analyzed) and the values 
of the parameters specified above (vq w A2wf for colli- 
sions with the incoming wall), B = 10~ 4 and rj is typ- 
ically 50%. This means that after three more collisions 
the particle leaves with an energy not enough to over- 
come the height of one single particle in the gravity field. 
For 7„ = 10 3 and the other parameters kept constant, 
B = 10 -5 and rj has been reduced to 5%, resulting in that 
the number of collisions needed for the particle to have 
its kinetic energy reduced to the same residual fraction, 
has increased roughly by an order of magnitude. On the 
other hand, given the weak dependence of Eq. (j^) on the 
velocity, one expects that the transitions shown in Fig. |^ 
will depend also weakly on the amplitude of the shaking 
velocity. The reduction of the inelasticity rj by an order 
of magnitude seems enough for particles to "climb" the 
walls and develop the characteristic surface rolls observed 
in numerical simulations. 



4 Discussion 

We have shown that the value of the normal damping co- 
efficient influences the convective pattern of horizontally 
shaken granular materials. By means of molecular dynam- 
ics simulations in two dimensions we can reproduce the 
pattern observed in real experiments, which corresponds 
to a situation of comparatively high damping, character- 
ized by inelasticity parameters rj larger than 5%. For lower 
damping, the upper layers of the material develop addi- 
tional surface rolls as has been reported previously. As 
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normal damping decreases, the lower rolls descend and fi- 
nally disappear completely at inelasticities of the order of 
1%. 
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